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Large-eddy  simulation  of  stratocumulus-topped 
atmospheric  boundary  layers  with  dynamic 
subgrid-scale  models 

By  Inane  Senocak 


1.  Motivation  and  objectives 

Earth’s  climate  and  its  geographical  variation  is  strongly  influenced  by  cloud  coverage. 
It  is  estimated  that  about  50%  of  the  earth  is  covered  by  clouds  at  any  given  time, 
providing  a shield  from  solar  radiation.  Radiative  energy  transfer  and  its  interaction 
with  clouds  play  an  important  role  in  the  thermal  structure  and  stratification  of  the 
atmosphere.  For  instance,  clouds  have  high  reflectivity  in  the  visible  wavelengths,  thus 
providing  relative  cooling  of  the  atmosphere.  They  also  absorb  strongly  in  the  infrared 
wavelengths,  resulting  in  heating  of  the  atmosphere  (Salby  1996). 

Condensation  is  the  major  physical  process  that  is  responsible  for  cloud  formation. 
Clouds  can  be  classified  into  four  broad  categories,  namely:  cumulus,  cirrus,  nimbus  and 
stratus  (Rogers  & Yau  1989).  Many  other  classifications  can  be  derived  from  combina- 
tions  of  these  four  broad  categories.  A comprehensive  description  can  be  found  in  Scorer 
& Wexler  (1967).  Among  various  types  of  clouds,  marine  stratocumulus  clouds  have  re- 
ceived increased  attention  because  of  the  important  role  they  play  in  the  global  radiation 
budget.  Marine  stratocumulus  clouds  cover  about  25%  of  the  Earth’s  ocean  at  any  in- 
stant. These  are  low-level  clouds  that  exist  below  1.5  km  with  several  hundred  meters  in 
thickness  and  they  rarely  produce  precipitation.  Their  horizontal  coverage  is  extensive 
and  more  homogeneous  than  other  type  of  clouds.  Their  appearance  is  grey  and  a wavy 
undersurface  is  typical  (Kantha  & Clayson  2000;  Heymsfield  1993;  Mason  1975). 

The  structure  of  the  marine  stratocumulus  cloud-topped  atmospheric  boundary  layer 
is  driven  by  both  cloud-top  radiative  cooling  and  positive  buoyancy  flux  from  the  surface 
that  maintains  the  atmospheric  boundary  layer  in  a well-mixed  turbulent  state.  Above 
the  cloud  layer,  negative  buoyancy  flux  has  a stabilizing  effect,  suppressing  the  turbulence 
there  (Kantha  & Clayson  2000).  The  entrainment  of  dry  air  from  above  the  cloud  layer 
induces  evaporative  cooling  and  entrained  air  parcel  can  sink  further  down  enhancing  the 
turbulent  mixing  within  the  clouds,  known  as  cloud  top  entrainment  instability  (Deardorff 
1980a).  Clearly,  the  structure  of  cloud- topped  atmospheric  boundary  is  more  complicated 
than  a cloud-free  atmospheric  boundary  layer  due  to  strong  interactions  among  cloud 
microphysics,  turbulent  motions  and  the  radiative  energy  transfer. 

In  large-eddy  simulation  (LES),  three-dimensional,  large  unsteady  flow  structures  are 
resolved  and  the  effects  of  the  unresolved  scales  are  modeled.  A filtering  operation  is 
applied  to  the  governing  equations  to  distinguish  between  the  resolved  scales  that  are 
computed  and  smaller  scales  that  are  modeled.  LES  has  been  widely  applied  to  simulate 
atmospheric  boundary  layers,  partly  because  of  the  difficulties  involved  in  observational 
studies  and  field  experiments  (Stevens  & Lenschow  2001).  An  extensively  studied  topic 
is  LES  of  cloud-free  convective  atmospheric  boundary  layers.  The  salient  features  of 
these  boundary  layers  have  been  compared  to  observations  and  well  documented  (Kantha 
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& Clayson  2000).  The  presence  of  clouds  complicates  the  problem  due  to  additional 
physics  and  reliable  numerical  simulations  of  cloud-topped  boundary  layer  is  still  an 
active  research  area.  In  the  following,  several  representative  studies  are  briefly  mentioned 
to  help  describe  the  current  state  of  the  knowledge. 

An  early  study  on  three-dimensional  modeling  of  cloud-topped  atmospheric  boundary 
layer  is  by  Deardorff  (19806).  He  studied  the  structure  of  turbulence  and  entrainment 
within  stratocumulus  layers  with  and  without  cloud-top  radiative  cooling  and  suggested 
that  simulations  need  high  resolution  at  the  inversion.  It  was  also  found  that  a func- 
tional dependence  on  Richardson  number  helps  correlate  the  entrainment  rate.  Deardorff 
(1980a)  has  defined  a criterion  for  the  cloud  top  entrainment  instability.  It  was  found 
that  entrainment  rate  increases  decisively  when  the  equivalent  potential  temperature  gra- 
dient across  the  cloud  top  drops  below  a critical  value.  It  was  also  shown  that  a strong 
instability  can  lead  to  stratocumulus  breakup  leading  to  a scattered  cumulus  layer.  Tag 
& Payne  (1987)  indicated  that  in  addition  to  DeardorfF’s  criterion,  the  vertical  motions 
should  exceed  a threshold  for  the  cloud  breakup  to  occur. 

Moeng  (1986)  studied  the  structure  of  a stratus- topped  boundary  layer  using  LES.  It 
was  found  that  the  vertical  component  of  the  turbulent  kinetic  energy  is  generated  by 
buoyancy  and  a portion  of  this  energy  is  redistributed  in  the  horizontal  directions  due 
to  pressure  effects.  It  was  also  shown  that  turbulence  is  generated  more  effectively  by 
surface  heating  than  cloud-top  cooling. 

Bohnert  (1993)  tested  the  dynamic  procedure  for  LES  of  cloud-topped  boundary  layer. 
Simple  parametrization  for  cloud  microphysics  and  radiation  were  adopted.  The  dynamic 
model  results  were  compared  to  the  results  of  SGS  model  that  were  optimized  in  an  ad- 
hoc  fashion.  The  results  were  found  to  be  comparable.  The  importance  of  SGS  modeling 
in  predicting  cloud  breakup  was  also  highlighted. 

Stevens  et  al.  (2000)  investigated  the  dependence  of  an  LES  model  on  mesh  resolution, 
numerical  schemes  and  SGS  model.  They  provided  simulations  of  varying  resolutions  for 
the  simulation  of  stratocumulus  topped  marine  boundary  layer.  Different  SGS  models 
and  advection  schemes  were  tested.  It  was  found  that  thickness  of  the  inversion  layer, 
depth  of  entraining  eddies  and  shape  of  the  vertical  velocity  spectra  is  influenced  by 
mesh  resolution.  Motions  at  the  inversion  were  found  to  be  underresolved  even  for  the 
finest  resolution.  The  entrainment  rate  was  found  to  depend  on  both  numerical  and  SGS 
dissipation. 

Duynkerke,  Zhang  & Jonker  (1995)  performed  an  observational  study  to  describe  the 
microphysical  and  turbulent  structure  of  stratocumulus  observed  during  the  Atlantic 
Stratocumulus  Transition  Experiment  (ASTEX).  The  turbulence  kinetic  energy  budget, 
velocity  and  temperature  variance,  and  vertical  fluxes  were  calculated.  The  entrainment 
was  found  to  be  very  efficient,  which  resulted  in  reduction  of  turbulent  kinetic  energy 
production  due  to  buoyancy.  It  was  also  shown  that  water  vapor  flux,  liquid  water  flux, 
and  drizzle  rate  have  the  same  magnitude. 

Stevens  et  al.  (1998)  presented  an  LES  study  of  the  ASTEX  case.  They  adopted  a 
drop-size  resolving  cloud  microphysics  model  that  enabled  them  to  perform  simulations 
with  and  without  drizzle.  It  was  found  that  inclusion  of  drizzle  in  modeling  leads  to 
sharp  decrease  in  entrainment  and  turbulent  kinetic  energy  generation  by  buoyancy.  The 
authors  have  hypothesized  that  shallow,  well-mixed  radiatively  driven  stratocumulus 
cannot  persist  in  the  presence  of  heavy  drizzle. 

Duynkerke  et  al.  (1999)  did  a comparison  of  actual  ASTEX  observations  with  com- 
putations obtained  from  LES  and  one-dimensional  single  column  models.  The  buoyancy 
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flux  obtained  from  LES  agrees  well  with  the  observations.  The  authors  concluded  that 
drizzle  has  small  influence  on  the  buoyancy  flux,  although  significant  uncertainty  exists 
in  its  parametrization. 

The  objective  of  the  present  study  is  to  evaluate  the  dynamic  procedure  in  LES  of 
stratocumulus  topped  atmospheric  boundary  layer  and  assess  the  relative  importance  of 
subgrid-scale  modeling,  cloud  microphysies  and  radiation  modeling  on  the  predictions. 
The  simulations  will  also  be  used  to  gain  insight  into  the  processes  leading  to  cloud  top 
entrainment  instability  and  cloud  breakup.  In  this  report  we  document  the  governing 
equations,  numerical  schemes  and  physical  models  that  are  employed  in  the  Goddard 
Cumulus  Ensemble  model  (GCEM3D).  We  also  present  the  subgrid-scale  dynamic  pro- 
cedures that  have  been  implemented  in  the  GCEM3D  code  for  the  purpose  of  the  present 
study. 


2.  Numerical  formulation 

The  numerical  model  used  in  this  study  to  simulate  cloud-topped  atmospheric  bound- 
ary layers  is  the  Goddard  Cumulus  Ensemble  model(GCEM3D).  Its  main  features  axe 
described  in  Tao  & Simpson  (1993),  Simpson  & Tao  (1993)  and  Tao  et  al.  (2003). 

2.1.  Governing  equations 

Acoustic  waves  are  part  of  the  solution  of  the  compressible  Navier-Stokes  equations  and 
very  small  time  steps  are  needed  to  resolve  them.  On  the  other  hand,  acoustic  waves  do 
not  impact  the  dynamics  of  thermal  convection  for  low  Mach  number  flows.  Therefore,  the 
governing  equations  of  motion  are  simplified  by  filtering  out  the  sound  waves  from  them. 
The  resulting  set  of  equations  is  known  as  the  anelastic  equations  (Ogura  & Phillips  1962). 
In  deriving  these  equations  the  basic  assumption  is  to  decompose  the  thermodynamic 
state  variables  into  a horizontally  averaged  base  quantity,  which  only  depends  on  altitude, 
and  a perturbation  quantity,  which  depends  on  both  time  and  space,  as  follows 

p(x,y,z,t)  = p0(z)  +p'{x,y,z,t), 
p(x,y,z,t)  = p0(z)  + p'(x,y,z,t), 

T(x,  y,  2,  t ) = Ta{z ) + T'(x,  y,  z,  t),  (2.1) 

where  p is  the  pressure,  p is  the  density  of  moist  air  and  T is  the  temperature.  The 
following  moist  equation  of  state  is  used  in  deriving  the  momentum  equations 

p = pRT{  1 + 0.61g„),  (2.2) 

where  R is  the  gas  constant  for  dry  air  and  qv  is  the  mixing  ratio  of  water  vapor.  The 
pressure  is  nondimensionalized  according  to  a reference  pressure  (p0)  value 

7r  = (P.)i?/cP(  (2.3) 

Po 

where  Cp  is  the  specific  heat  of  dry  air  at  constant  pressure.  The  potential  temperature 
is  defined  as 

7T 

The  anelastic  form  of  the  filtered  continuity  equation  reads  as  follows 

d(poUi)  n 


(2.4) 
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The  filtered  momentum  equations,  using  the  f-plane  approximation,  are  written  as 


diii  ^ 1 d(p0UiUj) 


dir'  dr, 


,9‘ 


—Cp90-^—  + — — I-  0.61^  — qi)  + FCOrioiis,  (2.6) 
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dt  ' p0  dxj  v dxi  dxj 

where  is  the  gravitational  acceleration.  The  Coriolis  term  that  appears  in  the  momen- 
tum equations  is  written  as 


Fcorioiis  — [2u)siruj)U2,  —2usin<j)ui,0], 


(2.7) 


where  to  is  the  angular  velocity  of  the  Earth  and  q>  gives  the  latitude,  r y represents  the 
subgrid  scale  stress  tensor  and  its  particular  form  is  described  in  section  2.3.  The  primed 
variables  represent  deviation  from  horizontally  averaged  quantities. 

The  equations  for  potential  temperature  0 and  the  water  vapor  mixing  ratio  qv  are 
written  as  follows 

ec  - er)  + %£(fa  + fg-rns-  mg) 

Op 

~hy^~{dice  sice ) "b  Qradi  (2-8) 


99  1 d(p06uj ) _ _aqf  Ly 

dt  p0  dxi  dxi  Cp 


dqv  1 d{p0qvUi)  dH 


dt 


+ 


Po 


dxi 


dxi 


-f  (c  ec  er)  + {d{ce  &ice)) 


(2.9) 


where  Lv,  Lf  and  Ls  are  the  latent  heats  of  condensation,  fusion  and  sublimation,  respec- 
tively. The  quantities  c,ec,er,  f,m,dice,SiCe  represent  the  rates  of  condensation,  evapo- 
ration of  cloud  droplets,  evaporation  of  rain  drops,  freezing  of  rain  drops,  melting  of  snow 
and  graupel/hail,  deposition  of  ice  particles  and  sublimation  of  ice  particles,  respectively. 
The  particular  forms  for  these  phase  change  rates  are  not  explicitly  formulated.  Instead, 
a saturation  adjustment  scheme  is  adopted  that  calculates  the  amount  of  phase  change 
rate  in  order  to  remove  any  supersaturated  vapor  and/or  subsaturation  of  cloud  water. 
The  saturation  adjustment  scheme  is  described  in  Soong  & Ogura  (1973)  and  Tao,  Simp- 
son &:  McCumber  (1989).  7 f is  the  subgrid  scale  flux  of  the  scalar,  which  is  explained  in 
section  2.3.  Qrad  is  the  source  due  to  radiative  heat  transfer.  It  is  described  in  section 
2.4. 


2.2.  Cloud  microphysics  model 

The  formulation  of  the  cloud  microphysical  processes  is  based  on  solving  scalar  transport 
equations  for  each  hydrometeor  species.  The  transport  equations  for  cloud  water  ( qc ), 
cloud  ice  (qice) , rain  (qr),  snow  (qs),  graupel/hail  (qg)  axe  written  as 


dqc  1 d(p0qcUi ) 

dt  Po  dxi 


dlic 

dxi 


+ (c  - ec)  + Tqc, 


(2.10) 


dqice  1 d{p0qiceiii ) 

dt  p0  dxi 

dqT  1 d{p0qrUi)  _ 1 d{p0qrUr ) 
ot  po  dxi  p0  dx3 

dq3  + 1 d{p0qsUj)  _ 1 d{p0qsUs) 
dt  po  dxi  po  dx3 


dxi 


~b  {dice 


~ Sice)  + Tqi 


Qice  5 


ft  Qr 

+ {ms  +mg-f3-fg-  er)  + TQr, 
'■  — + ( ds  — ss  — ms  + fs ) + Tqa, 


(2.11) 

(2.12) 

(2.13) 
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1 d(p0qgUi ) 1 d(p0qgUg ) _ ,,  _ , r \ , /0 

9f  p0  dr*  p0  dz3  024 

where  Fr,  Vs  and  Vg  are  the  fall  speeds  of  rain,  snow  and  graupel,  respectively.  Their 
values  are  obtained  from  parameterizations.  7 f is  the  subgrid-scale  flux  of  the  scalar. 
represents  the  microphysical  transfer  rates  between  the  hydrometeor  species.  A total  of  27 
different  processes  are  considered  and  the  reader  is  referred  to  Lin,  Farley  Orville  (1983) 
for  details  of  their  formulation.  To  illustrate  the  parameterizations  of  these  processes,  only 
the  microphysical  transfer  rate  of  rain  Tgr  is  described  in  this  report.  For  instance,  if  the 
temperature  is  above  0 °C,  then  the  production  term  for  rain  is  given  by  the  following 
equation 

Tgr  = Pracw  + Praut  + Psacw  + Pgacw  — Pgmlt  - Psmlt  + Prevp(1~  ^i).  (2.15) 

The  terms  on  the  right  hand  side  of  the  above  equation  are  the  accretion  of  cloud  water 
by  rain,  autoconversion  of  cloud  water  to  form  rain,  accretion  of  cloud  water  by  snow, 
accretion  of  cloud  water  by  graupel,  melting  of  graupel  to  form  rain,  melting  of  snow  to 
form  rain,  evaporation  of  rain,  respectively.  The  accretion  of  rain  by  cloud  water  Pracw 
is  written  as 

nErwnoraqcT(S  4-  6)  . . 

Pracw  — ’ (2.1b) 

where  Erw  is  the  collection  efficiency,  which  is  assumed  to  be  1,  nor  is  the  intercept 
parameter  of  the  rain  drop  size  distribution,  F is  the  gamma  function,  Ar  is  the  slope 
parameter  in  rain  size  distribution  and  a,  b are  empirical  constants.  Exponential  size 
distributions  are  assumed  for  rain,  snow  and  graupel/hail.  The  explicit  forms  of  other 
transformation  rates  are  documented  in  Lin  et  al.  (1983). 

2.3.  Subgrid-scale  turbulence  models 

Three  turbulence  closure  models  will  be  considered  for  comparison.  These  are  the  subgrid- 
scale  kinetic  energy  model  of  Klemp  & Wilhelmson  (1978), which  is  the  base  turbulence 
model  in  the  GCEM3D  code,  the  dynamic  Smagorinsky  model  of  Germano  et  al.  (1991) 
and  the  localized  dynamic  Smagorinsky  model  of  Piomelli  & Liu  (1995).  The  imple- 
mentation of  dynamic  turbulence  models  by  Kirkpatrick  (2002)  has  been  coupled  to  the 
GCEM3D  code.  In  the  following  sections  these  models  are  briefly  explained. 


2.3.1.  Subgrid-scale  kinetic  energy  model 

A transport  equation  for  subgrid-scale  turbulent  kinetic  energy  is  solved,  which  is  then 
used  to  specify  the  eddy  viscosity.  The  influence  of  buoyancy  on  the  turbulent  motions 
is  also  modeled.  The  equation  for  subgrid-scale  turbulent  kinetic  energy  is  written  as 


dE  , 1 d(poEui)  t/6'  , 

+ !_ - gw'{—  + 0.61g' 

at  p0  dxj  60 


dUi  d 


dE_ 

m Q_  ) ' 


dxi 


^£3/2,  (2.17) 


A = ( AxAyAz )1^3, 

Km  = CmAE1/2, 

Tij  = 3 ESij  - KmifaT.  + (2-18) 

where  Km  is  the  turbulent  eddy  viscosity,  A is  the  filter  width.  The  empirical  constants 
Ce  and  Cm  have  the  values  of  0.7  and  0.2, respectively.  The  subgrid-scale  scalar  fluxes 
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4>  TS  d(t> 


(2.19) 


where  Kh/Km  — 3 is  used.  The  buoyancy  flux  in  a saturated  area  is  computed  as 

dgt 
dz  ’ 


+ 0.61^  - q[)  = - AKh ^ 


where 


A = 


1 •,  , 1.61eLgv 


ft  i i ei-2^  ’ 

0 1 + CpJ?fr 


and  e = 0.622. 

In  an  unsaturated  area,  the  buoyancy  flux  is  computed  as  follows: 


^+0W,-9i)=-&(‘f  + 0.61^). 


(2.20) 


(2.21) 


(2.22) 


2.3.2.  Dynamic  Smagorinsky  model 

Smagorinsky  model  (Smagorinsky  1963)  is  commonly  used  in  LES  to  model  the  subgrid 
scale  stresses.  It  is  based  on  eddy  viscosity  assumption  and  can  be  written  as 

Tij  - ~6ijrkk  = —2CA2\S\Sij,  (2.23) 


where  C is  a dimensionless  parameter  and  |5|  = ^J2S~S^.  The  filtered  strain  rate  is 
defined  as 


_ 1 , dui  duj . 


(2.24) 


Germano  et  al.  (1991)  have  introduced  the  dynamic  procedure  that  offers  advantages 
over  the  original  Smagorinsky  model.  In  this  model  the  parameter  C is  computed  at  each 
time  step  from  the  information  already  available  in  the  resolved  velocity  field.  The  basic 
formalism  behind  the  dynamic  Smagorinsky  model  is  explained  in  the  following. 

The  dynamic  model  involves  filters  of  different  sizes.  In  addition  to  the  grid  filter,  a 
test  filter  is  introduced.  The  subgrid  scale  stress  tensor  based  on  the  grid  filter  and  the 
test  filter  are  written  respectively  as 


7"ij  — UilLj  UiUj  j 

Tij  — UiUj  UiUj,  (2.25) 

where  the  symbols  overbar  and  the  hat  represent  the  grid  and  test  filtering  operations, 
respectively.  Applying  the  test  filter  to  and  subtracting  it  from  Ty  yields  the  following 
identity  (Germano,  1992) 

Ltij  — 7) j Tij  — UiUj  UiUj . (2.26) 

The  significance  of  this  identity  is  that  it  can  be  computed  from  the  large  eddy  field. 
Germano  et  al.  (1991)  have  utilized  this  identity  to  dynamically  compute  the  coefficient 
of  the  Smagorinsky  model  as  follows. 


where 
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otij  = -2A  *SSi:j, 

Pij  = -2A  2SSi:j.  (2.28) 


For  atmospheric  boundary  layers,  where  the  horizontal  directions  are  assumed  to  ho- 
mogenous, the  filtering  operation  is  applied  only  in  the  horizontal  direction  and  C is 
assumed  to  be  independent  of  the  homogeneous  directions  and  taken  out  of  the  filtering 
operator.  Equation  (2.27)  is  rearranged  to  the  following  form 

Lij  ~ d2j  Lkk  — CMij,  (2.29) 

where 

Mij  = oiij  - Pij.  (2.30) 

Following  the  method  described  in  Lilly  (1992),  the  coefficient  C is  computed  so  as  to 
minimize  the  sum  of  the  squares  of  the  residuals  of  equation  (2.29).The  numerator  and 
the  denominator  are  averaged  over  the  horizontal  (x,y)-plane. 


C(z,t)  = 


Ljj)xy 

(MklMkl)  xy 


(2.31) 


Once  C is  calculated,  the  subgrid  scale  stress  tensor,  given  in  equation  (2.23)  is  computed. 
In  a similar  approach,  the  subgrid-scale  flux  of  a scalar  (f>  can  be  computed  dynamically 
Moin  et  al.  (1991).  If  we  consider  the  following  eddy  diffusivity  subgrid-scale  model 


Vt  d<f> 
Prt  dXi  ’ 


(2.32) 


where  the  kinematic  eddy  viscosity  is  computed  with  the  dynamic  Smagorinsky  model 
as  vt  = CA2|S|.  The  dynamic  procedure  can  also  be  applied  to  compute  the  turbulent 
Prandtl  number  Prt.  The  subgrid-scale  scalar  flux  based  on  the  test  filter  scale  is  written 
as 


(2.33) 


The  test  and  grid  scale  fluxes  are  related  to  each  other  by  the  following  identity 


Pi  = Gi  - 7i  = <j)Ui  - <j>Ui 


Prt  dxi 

The  above  equation  can  be  recast  into  the  following  form 


cAy|  dj  | cA2lgj  8$ 


Prt  dxi ' 


(2.34) 


Pi  = -Jp—Ri- 
Prt 


(2.35) 


Following  the  suggestion  of  Lilly  (1992),  a least  squares  procedure  is  applied  to  compute 
the  value  of  Prt.  Note  that  the  value  of  C is  determined  earlier. 


1 _ 1 ( PjRj)xy 

Prt  C (RmRm)xy 


(2.36) 


After  calculating  the  Prt,  equation  (2.31)  is  used  to  compute  the  subgrid  scale  flux  of  a 
scalar. 
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2.3.3.  Approximate  localized  dynamic  Smagorinsky  model 

The  dynamic  Smagorinsky  model  is  not  general  for  flows  with  no  homogeneous  di- 
rection, because  of  the  need  for  averaging  in  the  homogeneous  directions.  Ghosal  et  al. 
(1995)  have  addressed  the  mathematical  inconsistencies  and  proposed  a new  formulation 
for  the  dynamic  procedure  that  makes  it  applicable  to  arbitrary  inhomogeneous  flows. 
They  have  also  provided  formal  justifications  for  the  ad-hoc  procedures  that  have  been 
adopted  in  the  early  versions  of  the  dynamic  model. 

The  mathematical  inconsistency  in  previous  versions  comes  from  the  fact  that  space 
and  time  dependent  coefficient  C is  simply  taken  out  of  the  filtering  operation.  In  the 
variational  formulation  of  Ghosal  et  al.  (1995),  C is  kept  inside  and  an  integral  equation 
needs  to  be  solved  at  each  time  step  to  determine  C.  This  method  is  referred  to  as  the 
dynamic  localization  model. 

The  solution  of  an  integral  equation  is  costly.  Piomelli  & Liu  (1995)  have  followed  a 
simpler  approach  and  proposed  the  approximate  localized  dynamic  model.  In  the  follow- 
ing this  method  is  briefly  described. 

Equation  (2.27)  is  recast  in  the  following  form. 


-Caij  = L%  + C*/3y, 


(2.37) 


where 

Lij  — Lij  — -SijLkk.  (2.38) 

An  estimated  value  C*  is  assumed,  which  is  the  value  of  C from  the  previous  time  step. 
Along  with  this  approximation,  a least  squares  procedure  is  applied  to  compute  C as 
given  below 


C = -- 


(Lij  + C*Pij)oti 


(2.39) 

&7nn&mn 

The  same  iterative  idea  applies  to  equation  (2.31)  to  compute  the  Prt  of  the  subgrid- 
scale  flux  of  a scalar  quantity. 


2.4.  Radiation  model 

In  an  effort  to  describe  the  physical  models  in  GCEM3D  code  in  a single  document,  we 
briefly  summarize  the  basic  formalism  of  the  radiation  modeling. 

The  plane  parallel  assumption  is  typically  adopted  to  model  the  radiative  energy  trans- 
fer. The  convenience  of  the  assumption  comes  from  the  fact  that  the  properties  of  the 
atmosphere  vary  sharply  with  height  due  to  its  stratification,  hence,  the  medium  is  re- 
garded as  horizontally  homogeneous  (Salby  1996). 

The  intensity  of  a radiation  pencil,  traversing  a distance  ds  along  the  direction  of  its 
propagation,  changes  due  to  absorbtion,  scattering  and  emission,  which  can  be  described 
by  the  following  equation  (Liou  2002). 

dl\  = -K\pl\ds  + jxpds,  (2.40) 

where  kx  is  the  mass  extinction  cross  section  due  to  absorbtion  and  scattering  and  j\ 
is  the  source  function  due  to  emission  and  scattering.  In  plane-parallel  assumption,  it  is 
convenient  to  define  an  optical  thickness  as 

/OO 

kxdz,  dz  = pds,  p = cosd,  (2-41) 

where  z is  the  vertical  direction  and  6 is  the  angle  between  the  path  of  radiation  and  the 
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vertical,  which  is  specifically  called  the  zenith  angle.  The  basic  equation,  describing  the 
radiative  energy  transfer  in  plane-parallel  atmospheres  is  written  as 

pdI^ih-^2  = J(r;  p,  <j>)  - J(r;  p,  <j>).  (2.42) 

The  source  function  J is  defined  as 

J{t\p,4>)  = ^ J J I(t;  p',  <f>')P(p,  cj>\  p\  <f>')dp' d<f>'  + 

±FoP(p,ct>-,-p0,cf o)e~r^°  - (1  -a>)B[T(r)],  (2.43) 

47T 

where  P(p,  <j>\  p',  ft)  is  the  phase  function,  which  gives  the  angular  distribution  of  scat- 
tered energy  as  a function  of  direction,  <f  is  the  azimuthal  angle,  d>  is  the  single  scat- 
tering albedo,  B[T(t)\  is  the  Planck’s  function  representing  the  blackbody  emission  of 
the  medium,  and  F0  is  the  solar  irradiance  at  the  top  of  atmosphere.  Upon  solution  of 
equation  (2.43),  the  flux  density  F\  and  the  total  flux  density  F are  computed  based  on 
the  following  definitions 

f l r oo 

FI\t)  = 2-k  ll[{T,±p)pdp  / F\d\.  (2.44) 

Jo  Jo 


The  heating  rate  due  to  radiation  that  appears  in  the  potential  temperature  equation, 
Qrad,  is  then  computed  as 

1 dF(z ) 


Qrad  — 


p0Cp  dz 


(2.45) 


Equation  (2.43)  is  an  integrodifferential  equation  and  its  numerical  solution  is  quite 
involved  due  to  sharp  variation  of  the  atmospheric  properties  with  height.  The  compu- 
tation of  radiation  fluxes  involves  spectral,  vertical  and  directional  integrations.  Because 
the  absorbtion  coefficient  varies  sharply  with  wave  number,  the  spectral  integration  is  the 
most  CPU  intensive  part.  However,  the  major  difficulty  comes  from  the  dependence  of 
the  absorbtion  coefficient  on  pressure  and  temperature.  Hence,  effective  parametrization 
is  an  important  part  in  the  numerical  solution  of  equation  (2.42).  Different  approaches 
are  adopted  depending  on  the  nature  of  the  radiation  problem.  A detailed  account  of 
atmospheric  radiation  modeling  can  be  found  in  Liou  (2002). 

The  radiative  transfer  model  that  is  incorporated  into  the  GCEM  code  is  documented 
in  Chou  & Suarez  (1996a, 6).The  radiation  scheme  can  model  the  absorbtion  due  to  water 
vapor,  CO 2,  03,  and  O2,  and  scattering  by  clouds,  aerosols  and  molecules.  Fluxes  are 
integrated  almost  over  the  full  spectrum,  ranging  from  0.175  pm  to  lOjUm.  The  radiation 
field  is  divided  into  three  distinct  regions.  In  the  ultraviolet  and  photosynthetically  active 
regions,  the  spectrum  is  divided  into  8 bands,  in  which  single  03  absorption  coefficient 
and  Rayleigh  scattering  coefficient  is  used  for  each  band.  The  infrared  region  is  divided 
into  three  bands  and  the  k-distribution  method  is  applied  with  ten  absorption  coefficients 
in  each  band.  The  5-Eddington  approximation  is  used  to  compute  the  reflection  and 
transmission  of  a cloud  and  aerosol-laden  layer.  Two-stream  adding  method  is  then  used 
to  compute  the  fluxes. 


2.5.  Numerical  schemes 

In  GCEM3D  code,  the  governing  equations  are  discretized  on  a staggered  grid.  The  mo- 
mentum equations  are  solved  using  the  second  order  central  difference  scheme  for  the 
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spatial  derivatives  and  the  second  order  accurate  leap-frog  scheme  for  temporal  deriva- 
tives. A time  smoother  is  adopted  to  avoid  the  time  splitting  problem  associated  with 
the  leap-frog  scheme.  Forward  time  differencing  and  the  multi-dimensional  positive  def- 
inite advection  transport  algorithm  (MPDATA)  of  Smolarkiewicz  & Grabowski  (1990) 
are  used  to  solve  the  scalar  transport  equations.  A direct  solver  (FFT)  is  utilized  for  the 
solution  of  the  pressure  Poisson  equation.  The  GCEM3D  code  has  been  parallelized  to 
run  on  the  SGI  Origin  clusters  using  OpenMP  (Jin  et  al.  2003). 

Lateral  boundary  conditions  are  periodic  and  free  slip  boundary  conditions  are  imposed 
on  the  top  and  bottom  boundaries  for  all  the  variables  except  the  vertical  velocity  that 
vanishes  at  the  top  and  bottom  boundary. 


3.  Future  work 

The  sugrid-scale  turbulence  models  described  in  the  previous  section  will  be  tested 
with  progressively  refined  mesh  resolutions,  adopting  the  ASTEX  field  experiment  as 
the  test  case.  The  relative  influence  of  subgrid-scale  modeling,  cloud  microphysics  and 
radiation  modeling  on  the  predictions  will  be  investigated.  The  simulations  will  also  be 
used  to  gain  insight  into  the  processes  leading  to  the  cloud  top  entrainment  instability 
and  cloud  breakup.  Additionally,  the  dynamic  modeling  subroutine  will  be  parallelized 
to  speed  up  the  overall  computation. 

Preliminary  computations  produced  comparable  results  of  wind  speeds,  soundings  of 
cloud  water  and  potential  temperature  among  the  subgrid-scale  models  considered.  How- 
ever, the  computations  also  indicated  problems  with  the  top  boundary  resulting  in  spuri- 
ous cloud  formation  and  unstable  stratification.  A possible  cure  to  the  problem  might  be 
to  use  a Rayleigh  absorbing  layer,  which  is  commonly  used  in  atmospheric  simulations, 
in  the  proximity  of  the  top  boundary  to  damp  the  vertically  propagating  gravity  waves. 

In  this  report  we  have  briefly  documented  the  governing  equations  and  physical  mod- 
els that  are  employed  in  the  GCEM3D  code.  We  have  also  described  the  subgrid-scale 
dynamic  procedures  that  have  been  introduced  to  the  GCEM3D  code  in  this  study.  We 
are  in  the  process  of  double  checking  the  implementation  and  the  results.  The  findings 
of  the  present  study  will  be  published  in  the  open  literature. 
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